home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1996 July / C-C++ Users Group Library July 1996.iso / vol_200 / 212_01 / floatr.c < prev    next >
Encoding:
C/C++ Source or Header  |  1980-01-01  |  29.6 KB  |  1,326 lines
  1. /* FLOATR.C    VERS:- 01.00  DATE:- 09/26/86  TIME:- 09:36:35 PM */
  2. /*
  3. %CC1 $1.C 
  4. */
  5. /* 
  6. Description:
  7.  
  8. Library of floating point functions (FLOAT.C, by Bob Mathias)
  9.     with additions (FLOAT+44.C, by L.C. Calhoun)
  10.     with added trig functions (CTRIG.C, by L.C. Calhoun)
  11.     with added log functions (CLOGS.C, by L.C. Calhoun)
  12.     with additional function(s) (by J.A. Rupley)
  13.  
  14. FLOAT.C from BDS C distribution package, version 1.5
  15. FLOAT+44.C, CTRIG.C, AND CLOGS.C are from cug disk "functions III".
  16.  
  17. Note: the v1.5 FLOAT.C has an updated version of _spr;
  18. replaced v1.44 std functions by v1.5 std functions
  19.  
  20. By J.A. Rupley, Tucson, Arizona
  21. Coded for BDS C compiler, version 1.50a
  22. */
  23.  
  24.  
  25.     /* page eject */
  26.  
  27. /* summary of functions 
  28.  
  29. (from FLOAT, by Mathias)
  30. int fpcomp(op1,op2)            return = (1,0,-1) || (op1 >,=,< op2)
  31. char *fpadd(result,op1,op2)        result = op1 + op2
  32. char *fpsub(result,op1,op2)        result = op1 - op2
  33. char *fpmult(result,op1,op2)        result = op1*op2
  34. char *fpdiv(result,op1,op2)        result = op1/op2
  35. char *itof(result,int#)            result = float(int#)
  36. char *atof(result,string#)        result = float(string#)
  37. char *ftoa(result,fpno)            result = string(fpno)
  38. char *itoa(result,int#)            result = string(int#)
  39. void printf(format)            standard C print function
  40. void _spr(line,&format)            float formatting function
  41.  
  42. (from FLOAT+44, by Calhoun)
  43. char *fpmag(result,fpno)        result = abs(fpno)
  44. char *fpchs(result,fpno)        result = - fpno
  45. char *fpasg(result,fpno)        result = fpno
  46. int ftoit(fpno)                return = truncate then to integer
  47. int ftoir(fpno)                return = round then to integer
  48.  
  49. (from CLOGS, by Calhoun)
  50. int exprange(fpno)            service function
  51. char *pi(result)            result = float(pi)
  52. char *expe(result,fpno)            result = e^(fpno)
  53. char *exp10(result,fpno)        result = 10^(fpno)
  54. char *log10(result,&int,fpno)        result = log10(abs(fpno)) //
  55.                             *int = sgn(fpno)
  56.  
  57. (from CTRIG, by Calhoun)
  58. char *degtorad(result,fpno)        result = rad(fpno-degrees)
  59. char *radtodeg(result,fpno)        result = deg(fpno-radians)
  60. char *sinev(result,fpno)        service function
  61. char *sine(result,fpno)            result = sin(fpno-radians)
  62. char *cosine(result,fpno)        result = cos(fpno-radians)
  63. char *tangent(result,fpno)        result = tangent(fpno-radians)
  64. char *atanev(result,fpno)        service function
  65. char *arctan(result,fpno)        result = arctan-radians(fpno)
  66. char *arctan2(result,quadrant,opside,adjside)
  67.         result = arctan-radians(opside/adjside) // quadrant = 1 to 4
  68.     
  69. (additions, by Rupley)
  70. char *sqrt(result,&int,fpno)        result = sqrt(abs(fpno)) //
  71.                             *int = sgn(fpno)
  72. */
  73.  
  74.  
  75.     /* page eject */
  76.  
  77. /* FLOAT.C ***************/
  78.  
  79. /*
  80.     Floating point package support routines
  81.  
  82.     Note the "fp" library function, available in DEFF2.CRL,
  83.     is used extensively by all the floating point number
  84.     crunching functions.
  85.  
  86.     -------------------------------------------------------------
  87.     Usage: After compiling your program, link with this library
  88.     by typing:
  89.  
  90.     A>clink <your program files> -f float <cr>
  91.     -------------------------------------------------------------
  92.  
  93.     NEW FEATURE: a special "printf" function has been included
  94.              in this source file for use with floating point
  95.              operands, in addition to the normal types. The
  96.              printf presented here will take precedence over
  97.              the DEFF.CRL version when "float" is specified
  98.              on the CLINK command line at linkage time.
  99.              Note that the "fp" function, needed by most of
  100.              the functions in this file, resides in DEFF2.CRL
  101.              and will be automatically collected by CLINK.
  102.  
  103.     All functions here written by Bob Mathias, except printf and
  104.     _spr (written by Leor Zolman.)
  105. */
  106.  
  107. #include <bdscio.h>
  108.  
  109. #define NORM_CODE    0
  110. #define ADD_CODE    1
  111. #define SUB_CODE    2
  112. #define MULT_CODE    3
  113. #define DIV_CODE    4
  114. #define FTOA_CODE    5
  115.  
  116. fpcomp(op1, op2)
  117. char *op1, *op2;
  118. {
  119.     char work[5];
  120.     fpsub(work, op1, op2);
  121.     if (work[3] > 127)
  122.         return (-1);
  123.     if (work[0] + work[1] + work[2] + work[3])
  124.         return (1);
  125.     return (0);
  126. }
  127.  
  128. fpnorm(op1) char *op1;
  129. {
  130.     fp(NORM_CODE, op1, op1);
  131.     return (op1);
  132. }
  133.  
  134. fpadd(result, op1, op2)
  135. char *result, *op1, *op2;
  136. {
  137.     fp(ADD_CODE, result, op1, op2);
  138.     return (result);
  139. }
  140.  
  141. fpsub(result, op2, op1)
  142. char *result, *op1, *op2;
  143. {
  144.     fp(SUB_CODE, result, op1, op2);
  145.     return (result);
  146. }
  147.  
  148. fpmult(result, op1, op2)
  149. char *result, *op1, *op2;
  150. {
  151.     fp(MULT_CODE, result, op1, op2);
  152.     return (result);
  153. }
  154.  
  155. fpdiv(result, op1, op2)
  156. char *result, *op1, *op2;
  157. {
  158.     fp(DIV_CODE, result, op1, op2);
  159.     return (result);
  160. }
  161.  
  162. atof(fpno, s)
  163. char fpno[5], *s;
  164. {
  165.     char *fpnorm(), work[5], ZERO[5], FP_10[5];
  166.     int sign_boolean, power;
  167.  
  168.     initb(FP_10, "0,0,0,80,4");
  169.     setmem(fpno, 5, 0);
  170.     sign_boolean = power = 0;
  171.  
  172.     while (*s == ' ' || *s == '\t')
  173.         ++s;
  174.     if (*s == '-')
  175.     {
  176.         sign_boolean = 1;
  177.         ++s;
  178.     }
  179.     for (; isdigit(*s); ++s)
  180.     {
  181.         fpmult(fpno, fpno, FP_10);
  182.         work[0] = *s - '0';
  183.         work[1] = work[2] = work[3] = 0;
  184.         work[4] = 31;
  185.         fpadd(fpno, fpno, fpnorm(work));
  186.     }
  187.     if (*s == '.')
  188.     {
  189.         ++s;
  190.         for (; isdigit(*s); --power, ++s)
  191.         {
  192.             fpmult(fpno, fpno, FP_10);
  193.             work[0] = *s - '0';
  194.             work[1] = work[2] = work[3] = 0;
  195.             work[4] = 31;
  196.             fpadd(fpno, fpno, fpnorm(work));
  197.         }
  198.     }
  199.     if (toupper(*s) == 'E')
  200.     {
  201.         ++s;
  202.         power += atoi(s);
  203.     }
  204.     if (power > 0)
  205.         for (; power != 0; --power)
  206.             fpmult(fpno, fpno, FP_10);
  207.     else
  208.         if (power < 0)
  209.         for (; power != 0; ++power)
  210.             fpdiv(fpno, fpno, FP_10);
  211.     if (sign_boolean)
  212.     {
  213.         setmem(ZERO, 5, 0);
  214.         fpsub(fpno, ZERO, fpno);
  215.     }
  216.     return (fpno);
  217. }
  218. ftoa(result, op1)
  219. char *result, *op1;
  220. {
  221.     fp(FTOA_CODE, result, op1);
  222.     return (result);
  223. }
  224.  
  225. itof(op1, n)
  226. char *op1;
  227. int n;
  228. {
  229.     int *p;
  230.     p = op1;
  231.     p[0] = 0;
  232.     p[1] = n;
  233.     op1[4] = 15;
  234.     fp(NORM_CODE, op1, op1);
  235.     return op1;
  236. }
  237.  
  238. itoa(str, n)
  239. char *str;
  240. {
  241.     char *sptr;
  242.     sptr = str;
  243.     if (n < 0)
  244.     {
  245.         *sptr++ = '-';
  246.         n = -n;
  247.     }
  248.     _uspr(&sptr, n, 10);
  249.     *sptr = '\0';
  250.     return str;
  251. }
  252.  
  253. /*
  254.     This is the special formatting function, which supports the
  255.     "e" and "f" conversions as well as the normal "d", "s", etc.
  256.     When using "e" or "f" format, the corresponding argument in
  257.     the argument list should be a pointer to one of the five-byte
  258.     strings used as floating point numbers by the floating point
  259.     functions. Note that you don't need to ever use the "ftoa"
  260.     function when using this special printf/sprintf combination;
  261.     to achieve the same result as ftoa, a simple "%e" format
  262.     conversion will do the trick. "%f" is used to eliminate the
  263.     scientific notation and set the precision. The only [known]
  264.     difference between the "e" and "f" conversions as used here
  265.     and the ones described in the Kernighan & Ritchie book is that
  266.     ROUNDING does not take place in this version...e.g., printing
  267.     a floating point number which happens to equal exactly 3.999
  268.     using a "%5.2f" format conversion will produce " 3.99" instead
  269.     of " 4.00".
  270. */
  271.  
  272. printf(format)
  273. char *format;
  274. {
  275.     int putchar();
  276.     _spr(&format, &putchar);        /* use "_spr" to form the output */
  277. }
  278.  
  279.  
  280. _spr(fmt, putcf, arg1)
  281. int (*putcf)();
  282. char **fmt;
  283. {
  284.     char _uspr(), c, base, *sptr, *format;
  285.     char wbuf[MAXLINE], *wptr, pf, ljflag, zfflag;
  286.     int width, precision, exp, *args;
  287.  
  288.     format = *fmt++;        /* fmt first points to the format string */
  289.     args = fmt;        /* now fmt points to the first arg value */
  290.     while (c = *format++)
  291.         if (c == '%')
  292.         {
  293.             wptr = wbuf;
  294.             precision = 6;
  295.             ljflag = pf = zfflag = 0;
  296.  
  297.             if (*format == '-')
  298.             {
  299.                 format++;
  300.                 ljflag++;
  301.             }
  302.  
  303.             if (*format == '0')
  304.                 zfflag++;        /* test for zero fill */
  305.  
  306.             width = isdigit(*format) ? _gv2(&format) : 0;
  307.  
  308.             if ((c = *format++) == '.')
  309.             {
  310.                 precision = _gv2(&format);
  311.                 pf++;
  312.                 c = *format++;
  313.             }
  314.  
  315.             switch (toupper(c))
  316.             {
  317.             case 'E' :
  318.                 if (precision > 7)
  319.                     precision = 7;
  320.                 ftoa(wbuf, *args++);
  321.                 strcpy(wbuf + precision + 3, wbuf + 10);
  322.                 width -= strlen(wbuf);
  323.                 goto pad2;
  324.  
  325.             case 'F' :
  326.                 ftoa(&wbuf[60], *args++);
  327.                 sptr = &wbuf[60];
  328.                 while (*sptr++ != 'E')
  329.                     ;
  330.                 exp = atoi(sptr);
  331.                 sptr = &wbuf[60];
  332.                 if (*sptr == ' ')
  333.                     sptr++;
  334.                 if (*sptr == '-')
  335.                 {
  336.                     *wptr++ = '-';
  337.                     sptr++;
  338.                     width--;
  339.                 }
  340.                 sptr += 2;
  341.  
  342.                 if (exp < 1)
  343.                 {
  344.                     *wptr++ = '0';
  345.                     width--;
  346.                 }
  347.  
  348.                 pf = 7;
  349.                 while (exp > 0 && pf)
  350.                 {
  351.                     *wptr++ = *sptr++;
  352.                     pf--;
  353.                     exp--;
  354.                     width--;
  355.                 }
  356.  
  357.                 while (exp > 0)
  358.                 {
  359.                     *wptr++ = '0';
  360.                     exp--;
  361.                     width--;
  362.                 }
  363.  
  364.                 *wptr++ = '.';
  365.                 width--;
  366.  
  367.                 while (exp < 0 && precision)
  368.                 {
  369.                     *wptr++ = '0';
  370.                     exp++;
  371.                     precision--;
  372.                     width--;
  373.                 }
  374.  
  375.                 while (precision && pf)
  376.                 {
  377.                     *wptr++ = *sptr++;
  378.                     pf--;
  379.                     precision--;
  380.                     width--;
  381.                 }
  382.  
  383.                 while (precision > 0)
  384.                 {
  385.                     *wptr++ = '0';
  386.                     precision--;
  387.                     width--;
  388.                 }
  389.  
  390.                 goto pad;
  391.  
  392.  
  393.             case 'D' :
  394.                 if (*args < 0)
  395.                 {
  396.                     *wptr++ = '-';
  397.                     *args = -*args;
  398.                     width--;
  399.                 }
  400.             case 'U' :
  401.                 base = 10;
  402.                 goto val;
  403.  
  404.             case 'X' :
  405.                 base = 16;
  406.                 goto val;
  407.  
  408.             case 'O' :
  409.                 base = 8;
  410.  
  411. val :
  412.                 width -= _uspr(&wptr, *args++, base);
  413.                 goto pad;
  414.  
  415.             case 'C' :
  416.                 *wptr++ = *args++;
  417.                 width--;
  418.                 goto pad;
  419.  
  420.             case 'S' :
  421.                 if (!pf)
  422.                     precision = 200;
  423.                 sptr = *args++;
  424.                 while (*sptr && precision)
  425.                 {
  426.                     *wptr++ = *sptr++;
  427.                     precision--;
  428.                     width--;
  429.                 }
  430.  
  431. pad :
  432.                 *wptr = '\0';
  433. pad2 :
  434.                 wptr = wbuf;
  435.                 if (!ljflag)
  436.                     while (width-- > 0)
  437.                         if ((*putcf)(zfflag ? '0' : ' ', arg1) == ERROR)
  438.                             return ERROR;
  439.  
  440.                 while (*wptr)
  441.                     if ((*putcf)(*wptr++, arg1) == ERROR)
  442.                         return ERROR;
  443.  
  444.                 if (ljflag)
  445.                     while (width-- > 0)
  446.                         if ((*putcf)(' ', arg1) == ERROR)
  447.                             return ERROR;
  448.                 break;
  449.  
  450.             default :
  451.                 if ((*putcf)(c, arg1) == ERROR)
  452.                     return ERROR;
  453.             }
  454.         }
  455.     else
  456.         if ((*putcf)(c, arg1) == ERROR)
  457.         return ERROR;
  458. }
  459.  
  460.  
  461.     /* page eject */
  462.  
  463. /* FLOAT+44.C ************/
  464. /*
  465.  
  466.     New Functions Added    fpmag converts to floating magnitude
  467.                   fpchs changes sign of floating point no
  468.                   fpasg provides assignment of fl pt no
  469.                   ftoit converts fl pt no to trunc. int.
  470.                   ftoir converts fl pt no to rounded int.
  471.  
  472.                 written by L. C. Calhoun
  473.     Second Revision by L. C. Calhoun  
  474.                 Modify the _spr to use the z option
  475.                   as in STDLIB2
  476.                 Modify the program set to utilize the
  477.                   V 1.44 zero insert string variable
  478.  
  479. */
  480.  
  481. #define EXPON_SIGN    0x80   /* break point for exponent sign */
  482.  
  483. /*  NEW FUNCTIONS ADDED */
  484. /*  new function to convert floating to truncated integer */
  485. ftoit(fpno)
  486. char *fpno;
  487. {
  488.     char *ftoa(), temp[20], temp2[20], *sptr, *wptr;
  489.     int atoi(), exp, pf;
  490.     wptr = &temp2;
  491.     ftoa(temp, fpno);
  492.     sptr = &temp;
  493.     while (*sptr++ != 'E')
  494.         ;
  495.     exp = atoi(sptr);
  496.     sptr = &temp;
  497.     if (*sptr == ' ')
  498.         sptr++;
  499.     if (*sptr == '-')
  500.     {
  501.         *wptr++ = '-';
  502.         sptr++;
  503.     }
  504.     sptr += 2;
  505.  
  506.     if (exp < 1)
  507.         *wptr++ = '0';
  508.  
  509.     pf = 7;
  510.     while (exp > 0 && pf)
  511.     {
  512.         *wptr++ = *sptr++;
  513.         pf--;
  514.         exp--;
  515.     }
  516.     while (exp > 0)
  517.     {
  518.         *wptr++ = '0';
  519.         exp--;
  520.     }
  521.     *wptr++ = '.';
  522.     /* foregoing lifted from _spr with F format  */
  523.     return (atoi(temp2));
  524. }
  525.  
  526. /* new function to round, then convert to integer */
  527. ftoir(fpno)
  528. char *fpno;
  529. {
  530.     char *fpsub(), *fpadd();
  531.     char *atof(), rnd[5], res[5];
  532.     int ftoit();
  533.     atof(rnd, "0.5");
  534.     if (fpno[3] > 127)
  535.         fpsub(res, fpno, rnd);
  536.     else
  537.         fpadd(res, fpno, rnd);
  538.     return (ftoit(res));
  539. }
  540.  
  541. /* new function to produce the unsigned magnitude of a fp no. */
  542. char *fpmag(result, fpno)
  543. char *result, *fpno;
  544. {
  545.     char *fpchs();
  546.     if (fpno[3] < 128)
  547.         return (fpasg(result, fpno));
  548.     else
  549.         return (fpchs(result, fpno));
  550. }
  551.  
  552. /* new function to change sign of floating point number */
  553. char *fpchs(result, fpno)
  554. char *result, *fpno;
  555. {
  556.     char *fpsub(), *ZERO;
  557.     ZERO = "\0\0\0\0\0";
  558.     return (fpsub(result, ZERO, fpno));
  559. }
  560.  
  561. /* a new function to assign <in essence> the value of one
  562.    floating point number to another  */
  563. char *fpasg(result, fpno)
  564. char *result, *fpno;
  565. {
  566.     int index;
  567.     for (index = 0; index <= 4; index++)
  568.         result[index] = fpno[index];
  569.     return (result);
  570. }
  571.  
  572.  
  573.     /* page eject */
  574.  
  575. /*  CLOGS.C       *****
  576. A group of programs in C, using the BDS-C floating point package
  577. as modified by LCC (FLOAT+44.*) and depending on the ability to
  578. insert null characters in a string available in BDS-C V 1.44.
  579. Four functions are handled:
  580.     log10, exp10, expe, pi
  581. In addition, there is a service function exprange() which returns
  582. a false (00) if the exponent of a floating point variable is reaching
  583. near the end of the usable range.
  584. These are discussed in detail in CLOGS.DOC
  585.  
  586. L. C. Calhoun
  587. 257 South Broadway
  588. Lebanon, Ohio 45036  513 day 433 7510 nite 932 4541
  589.  
  590.   */
  591.  
  592.  
  593. char *pi(result)
  594.  
  595. /* result is a standard character array char result[5]; in calling
  596. program as used for floating point variables.  The return is a pointer
  597. to result, and the value of pi stored in the result array in floating
  598. point */
  599.  
  600. char *result;
  601. {
  602.     char *piconst, *fpasg();
  603.  
  604.     piconst = "\171\356\207\144\2";
  605.     return (fpasg(result, piconst));
  606. }
  607.  
  608. char *expe(result, xin)
  609.  
  610. /* computes the base of the natural log "e" raised to the "x'th"
  611.     power.  Checks are made for out of range values and result is
  612.     defaulted to 0, 1, or a large number as appropriate.  There are
  613.     no error flags.  The arguments are floating point character
  614.     arrays char result[5], x[5]; in calling program.  Return is
  615.     a pointer to result, and "e to the x" stored in result.
  616.  */
  617.  
  618. char *result, *xin;
  619. {
  620.     char *zero, *one, *large, *coef[7], *eghty6, *meghty6;
  621.     char intres[5], xint[5], x[5];
  622.     char *fpmult(), *fpasg(), *fpdiv(), *fpchs();
  623.     int signx, index, bigexp, smallexp, zeroin;
  624.     int exprange();
  625.  
  626.     bigexp = smallexp = zeroin = 0;
  627.  
  628.     zero = "\0\0\0\0\0";
  629.     one = "\0\0\0\100\1";
  630.     large = "\0\0\0\100\175";
  631.     eghty6 = "\0\0\0\126\7";
  632.     meghty6 = "\0\0\0\252\7";
  633.  
  634.     coef[0] = "\0\0\0\100\1";
  635.     coef[1] = "\140\326\377\177\376";
  636.     coef[2] = "\130\373\3\100\374";
  637.     coef[3] = "\200\1\352\124\370";
  638.     coef[4] = "\351\253\362\131\364";
  639.     coef[5] = "\21\213\32\133\357";
  640.     coef[6] = "\371\330\260\134\354";
  641.  
  642.     /* preserve input datum */
  643.     fpasg(x, xin);
  644.  
  645.     /* check for sign */
  646.     if (x[3] < 128)        /* positive */
  647.     {
  648.         signx = 1;
  649.     }
  650.     else
  651.         {
  652.         signx = 0;
  653.         fpchs(x, x);
  654.     }
  655.  
  656.     /* check for zero and out of range of fp var */
  657.  
  658.     /* check for zero and very small numbers */
  659.     if (((x[4] > 127) && (x[4] < 226)) || ((x[4] == 0) &&
  660.         (x[3] == 0) && (x[2] == 0) && (x[1] == 0) && (x[0] == 0)))
  661.         zeroin = 1;
  662.  
  663.     /* check for very small exponent */
  664.     if (fpcomp(xin, meghty6) < 0)
  665.         smallexp = 1;
  666.  
  667.     /*  check for very large exponent */
  668.     if (fpcomp(x, eghty6) > 0)
  669.         bigexp = 1;
  670.  
  671.     /*  now if small number or zero, result is one */
  672.     /*  now if big number and positive, result is large number */
  673.     /*  now if big number and negative, result is zero */
  674.  
  675.     if (zeroin)
  676.         return (fpasg(result, one));
  677.     if (smallexp)
  678.         return (fpasg(result, zero));
  679.     if (bigexp)
  680.         return (fpasg(result, large));
  681.  
  682.     /*  all exceptions taken care of, so evaluate rest  */
  683.     fpasg(result, one);
  684.     fpasg(xint, x);
  685.     index = 1;
  686.     while ((index <= 6) && exprange(xint))
  687.     {
  688.         fpmult(intres, coef[index], xint);
  689.         fpadd(result, result, intres);
  690.         fpmult(xint, xint, x);
  691.         index++;
  692.     }
  693.     /* now do the square square */
  694.     fpmult(result, result, result);
  695.     fpmult(result, result, result);
  696.  
  697.     /* now treat sign appropriately */
  698.     if (signx)
  699.         return (result);
  700.     else
  701.         {
  702.         fpdiv(result, one, result);
  703.         return (result);
  704.     }
  705. }
  706.  
  707. char *exp10(result, xin)
  708.  
  709. /* similar to expe, except result returned is 10 raised to the x
  710. power    the antilogarithm to base 10  */
  711.  
  712. char *result, *xin;
  713. {
  714.     char *zero, *ten, *one, *large, *thty8;
  715.     char xint[5], *coef[7], intres[5], tenfac[5], x[5];
  716.     int index, bigexp, smallexp, signx, tenpower;
  717.     int exprange();
  718.  
  719.     bigexp = smallexp = 0;
  720.  
  721.     zero = "\0\0\0\0\0";
  722.     one = "\0\0\0\100\1";
  723.     large = "\0\0\0\100\175";
  724.     ten = "\0\0\0\120\4";
  725.     thty8 = "\0\0\0\114\6";
  726.  
  727.     coef[1] = "\65\264\256\111\1";
  728.     coef[2] = "\0\14\330\124\0";
  729.     coef[3] = "\0\46\354\100\377";
  730.     coef[4] = "\24\140\107\115\375";
  731.     coef[5] = "\242\304\361\155\372";
  732.     coef[6] = "\361\143\246\134\371";
  733.  
  734.     /* preserve input datum */
  735.     fpasg(x, xin);
  736.  
  737.     /* check for sign */
  738.     if (x[3] < 128)        /* positive */
  739.         signx = 1;
  740.     else        /* negative */
  741.     {
  742.         signx = 0;
  743.         fpchs(x, x);
  744.     }
  745.  
  746.     /* check for very small or large numbers, check by exponent size */
  747.     /* check for zero or small */
  748.     if (((x[4] > 127) && (x[4] < 226)) || ((x[4] == 0) &&
  749.         (x[3] == 0) && (x[2] == 0) && (x[1] == 0) && (x[0] == 0)))
  750.         smallexp = 1;
  751.     /* check for big number */
  752.     if (fpcomp(x, thty8) > 0)
  753.         bigexp = 1;
  754.  
  755.     /* if value is small or zero, return 1 as with expe */
  756.     /* if value is large and positive, return a large number */
  757.     /* if value is large and negative, return zero */
  758.  
  759.     if (smallexp)
  760.         return (fpasg(result, one));
  761.  
  762.     if (bigexp && signx)
  763.         return(fpasg(result, large));
  764.  
  765.     if (bigexp && !signx)
  766.         return(fpasg(result, zero));
  767.  
  768.     /* now reduce range of x to between zero and one */
  769.  
  770.     tenpower = ftoit(x);
  771.     itof(tenfac, tenpower);
  772.     fpsub(x, x, tenfac);
  773.     fpasg(tenfac, one);
  774.     while (tenpower)
  775.     {
  776.         fpmult(tenfac, tenfac, ten);
  777.         tenpower--;
  778.     }
  779.     /* now evaluate series  */
  780.     fpasg(result, one);
  781.     fpasg(xint, x);
  782.     index = 1;
  783.  
  784.     while ((index <= 6) && exprange(xint))
  785.     {
  786.         fpmult(intres, coef[index], xint);
  787.         fpadd(result, result, intres);
  788.         fpmult(xint, xint, x);
  789.         index += 1;
  790.     }
  791.  
  792.     /* now square result (note error in referenced article) */
  793.     fpmult(result, result, result);
  794.  
  795.     /* now check sign and make proper return */
  796.     fpmult(result, result, tenfac);        /* scale back up */
  797.  
  798.     if (signx)
  799.         return (result);
  800.  
  801.     else
  802.         return (fpdiv(result, one, result));
  803. }
  804.  
  805.  
  806. char *log10(result, sign, xin)
  807.  
  808. /* computes briggsian logarithm of x which is a char[5]
  809.     floating point number.  Return is logarithm in result[5],
  810.     and sign pointed to by sign.  The logarithm is taken
  811.     of the magnitude, and sign information preserved
  812.     as required by sign.
  813.  */
  814. char *result, *xin;
  815. int *sign;
  816. {
  817.     char *zero, *ten, *one, *large;
  818.     char *sqrtten, x[5];
  819.     char gamma[5], gamnum[5], gamden[5], *coef[5];
  820.     char *half;
  821.     char intres[5], gamint[5];
  822.     int tenpower;
  823.     int index, bigexp, smallexp, signx;
  824.     int exprange();
  825.  
  826.     bigexp = smallexp = 0;
  827.  
  828.     zero = "\0\0\0\0\0";
  829.     one = "\0\0\0\100\1";
  830.     large = "\0\0\0\114\6";
  831.     ten = "\0\0\0\120\4";
  832.     half = "\0\0\0\100\0";
  833.     sqrtten = "\304\145\61\145\2";
  834.  
  835.     coef[0] = "\362\6\56\157\0";
  836.     coef[1] = "\30\346\21\112\377";
  837.     coef[2] = "\100\55\344\132\376";
  838.     coef[3] = "\106\73\244\140\375";
  839.     coef[4] = "\174\5\367\141\376";
  840.  
  841.     /*  preserve input variable */
  842.     fpasg(x, xin);
  843.  
  844.     /* check for sign */
  845.     if (x[3] < 128)        /* positive */
  846.         signx = 1;
  847.     else        /* negative */
  848.     {
  849.         signx = -1;
  850.         fpchs(x, x);
  851.     }
  852.  
  853.     /* check for very small or large numbers, check by exponent size */
  854.     /* check for zero or small */
  855.     if (((x[4] > 127) && (x[4] < 209)) || ((x[4] == 0) &&
  856.         (x[3] == 0) && (x[2] == 0) && (x[1] == 0) && (x[0] == 0)))
  857.         smallexp = 1;
  858.     /* check for big number */
  859.     if ((x[4] > 47) && (x[4] < 128))
  860.         bigexp = 1;
  861.  
  862.     /* if very small, return - a large number
  863.         if very large, return + a large number  */
  864.  
  865.     if (smallexp)
  866.     {
  867.         *sign = signx;
  868.         return (fpchs(result, large));
  869.     }
  870.     if (bigexp)
  871.     {
  872.         *sign = signx;
  873.         return (fpasg(result, large));
  874.     }
  875.     /*  now bring into range 1 <= x < 10 */
  876.     tenpower = 0;
  877.     while (fpcomp(x, ten) >= 0)
  878.     {
  879.         tenpower++;
  880.         fpdiv(x, x, ten);
  881.     }
  882.     while (fpcomp(x, one) < 0)
  883.     {
  884.         tenpower--;
  885.         fpmult(x, x, ten);
  886.     }
  887.  
  888.     /* now if exactly one, no need to evaluate */
  889.     fpsub(gamnum, x, one);
  890.     if (((gamnum[4] > 127) && (gamnum[4] < 209)) || ((gamnum[0] == 0) &&
  891.         (gamnum[1] == 0) && (gamnum[2] == 0) && (gamnum[3] == 0)))
  892.     {
  893.         *sign = signx;
  894.         itof(result, tenpower);
  895.         return (result);
  896.     }
  897.     /* now compute gamma  for series  */
  898.  
  899.     fpsub(gamnum, x, sqrtten);
  900.     /* now check for size of numerator */
  901.     if (((gamnum[4] > 127) && (gamnum[4] < 209)) || ((gamnum[0] == 0) &&
  902.         (gamnum[1] == 0) && (gamnum[2] == 0) && (gamnum[3] == 0)))
  903.     {
  904.         itof(result, tenpower);
  905.         fpadd(result, result, half);
  906.         *sign = signx;
  907.         return (result);
  908.     }
  909.     fpadd(gamden, x, sqrtten);
  910.     fpdiv(gamma, gamnum, gamden);
  911.  
  912.     /* now set up for series (use gamnum as gamma squared) */
  913.     fpmult(gamnum, gamma, gamma);
  914.     fpasg(gamint, gamma);
  915.     index = 0;
  916.     fpasg(result, half);
  917.  
  918.     /* now do series evaluation */
  919.     while ((index <= 4) && exprange(gamint))
  920.     {
  921.         fpmult(intres, coef[index], gamint);
  922.         fpadd(result, result, intres);
  923.         fpmult(gamint, gamint, gamnum);
  924.         index++;
  925.     }
  926.  
  927.     /* now do correction for range reduction */
  928.     if (tenpower != 0)
  929.     {
  930.         itof(intres, tenpower);
  931.         fpadd(result, result, intres);
  932.     }
  933.  
  934.     /* now clean up and return */
  935.  
  936.     *sign = signx;
  937.     return (result);
  938. }
  939.  
  940. int exprange(x)
  941.  
  942. /* The input argument is a floating point function from BDS C which
  943. consists of an array of 5 character data.  The function returns
  944. a 1 if the exponent is in the range of - 47 to + 47.  Outside this
  945. range a value of 0 (false) is returned.  This is a range of ten to
  946. plus or minus 14 power */
  947. char *x;
  948.  
  949. {
  950.     if (((x[4] < 128) && (x[4] > 47)) || ((x[4] > 127) && (x[4] < 209)))
  951.         return (0);
  952.     else
  953.         return (1);
  954. }
  955.  
  956.  
  957.     /* page eject */
  958.  
  959. /*  **** CTRIG.C *****
  960. A group of programs in C, using the BDS-C floating point package,
  961. as modified by L. C. Calhoun called FLOATXT, which compute some
  962. commonly used transcendental functions - to wit*
  963.     sine, cosine, tangent and arctangent
  964.     convert degrees to radians, convert radians to degrees
  965. These functions are discussed in detail in CTRIG.DOC
  966.  
  967. L. C. Calhoun
  968. 257 South Broadway
  969. Lebanon, Ohio 45036   513 932 4541/433 7510
  970.  
  971.  *** revision 13 July 1981 to 1) refine precision of constants,
  972.  especially pi related values  2) to utilize string constants
  973.  as pseudo static numeric constants.  THIS WILL ONLY WORK WITH
  974.  BDS C V 1.44 (and hopefully later).
  975.  3) to add ARCTAN2() function to identify quadrant
  976.  
  977.   */
  978.  
  979. /* simple ones first converting degrees - radians */
  980.  
  981. char *degtorad(rad, deg)        /*obvious arguments in 5 char fp */
  982. char *rad, *deg;
  983. {
  984.     char *fpmult(), *radindeg;
  985.     radindeg = "\71\36\175\107\373";
  986.     fpmult(rad, deg, radindeg);
  987.     return (rad);
  988. }
  989.  
  990. char *radtodeg(deg, rad)        /* 5 char fp arguments */
  991. char *deg, *rad;
  992. {
  993.     char *fpmult(), *deginrad;
  994.     deginrad = "\12\162\227\162\6";
  995.     fpmult(deg, rad, deginrad);
  996.     return (deg);
  997. }
  998.  
  999. /* service function sinev which evaluates when range of angle
  1000. reduced to plus or minus pi/2 (90 deg) */
  1001.  
  1002. char *sinev(result, angle)
  1003.  
  1004. char *result, angle[5];
  1005. {
  1006.     char *fpmult(), x[5], xsq[5];
  1007.     char *coef[5], termreslt[5];
  1008.     char *fpadd(), *fpasg();
  1009.     int index;
  1010.  
  1011.     /*  use the exponent part of the floating point
  1012.             to check for threat of underflow  use small
  1013.             angle approximation if appropriate  */
  1014.     if ((angle[4] > 128) && (angle[4] < 226))
  1015.     {
  1016.         fpasg(result, angle);
  1017.         return (result);
  1018.     }
  1019.     /* solution to fpmult underflow problem */
  1020.  
  1021.     /* series coef are 1., -.1666666, .008333026, -.0001980742,
  1022.         .000002601887  determined from coefset program */
  1023.     coef[0] = "\0\0\0\100\1";
  1024.     coef[1] = "\157\252\252\252\376";
  1025.     coef[2] = "\271\242\103\104\372";
  1026.     coef[3] = "\320\352\46\230\364";
  1027.     coef[4] = "\246\15\116\127\356";
  1028.     fpasg(x, angle);
  1029.     fpmult(xsq, x, x);
  1030.     setmem(result, 5, 0);
  1031.     /* to this point the coef have been initialized, the angle
  1032.         copied to x, x squared computed, and the result initialized */
  1033.  
  1034.     /* now to do the polynomial approximation */
  1035.     index = 0;
  1036.     while ((index <= 4) && ((x[4] > 194) || (x[4] < 64)))
  1037.         /* use index for loop, and exponent of x to avoid underflow
  1038.                problems */
  1039.     {
  1040.         fpmult(termreslt, coef[index], x);
  1041.         fpadd(result, result, termreslt);
  1042.         index++;
  1043.         fpmult(x, x, xsq);
  1044.     }
  1045.     return (result);
  1046. }
  1047.  
  1048. /* here is sine(result,angle) with angle in radians */
  1049.  
  1050. char *sine(result, angle)
  1051.  
  1052. char *result, *angle;
  1053. {
  1054.     char *fpmult(), *twopi, *halfpi;
  1055.     char *mtwopi, *mhalfpi, *fpasg(), *fpchs();
  1056.     char *pi, *sinev(), *fpadd();
  1057.     char y[5], *fpsub();
  1058.     int fpcomp(), compar;
  1059.     int signsine;
  1060.     /* some initialization required here */
  1061.     signsine = 1;
  1062.     twopi = "\171\356\207\144\3";
  1063.     halfpi = "\171\356\207\144\1";
  1064.     pi = "\171\356\207\144\2";
  1065.     mtwopi = "\207\21\170\233\3";
  1066.     mhalfpi = "\207\21\170\233\1";
  1067.     fpasg(y, angle);
  1068.     while (fpcomp(y, twopi) >= 0)
  1069.     {
  1070.         fpsub(y, y, twopi);
  1071.     }
  1072.     while (fpcomp(y, mtwopi) <= 0)
  1073.     {
  1074.         fpadd(y, y, twopi);
  1075.     }
  1076.     if (fpcomp(y, halfpi) > 0)
  1077.     {
  1078.         signsine *= -1;
  1079.         fpsub(y, y, pi);
  1080.     }
  1081.     if (fpcomp(y, mhalfpi) < 0)
  1082.     {
  1083.         signsine *= -1;
  1084.         fpadd(y, y, pi);
  1085.     }
  1086.     sinev(result, y);
  1087.     if (signsine > 0)
  1088.         return(result);
  1089.     /* minus so need to change sign */
  1090.     else
  1091.         return (fpchs(result, result));
  1092. }
  1093.  
  1094. /* cosine(result,angle) with angle in radians  - uses sine */
  1095.  
  1096. char *cosine(result, angle)
  1097.  
  1098. char *result, *angle;
  1099. {
  1100.     char *sine(), *fpsub(), *halfpi, y[5];
  1101.     halfpi = "\171\356\207\144\1";
  1102.     fpsub(y, halfpi, angle);
  1103.     sine(result, y);
  1104.     return (result);
  1105. }
  1106.  
  1107. /* tangent(result,angle) with angle in radians, returns very 
  1108. large number for divide by zero condition */
  1109.  
  1110. char *tangent(result, angle)
  1111. char *result, angle[5];
  1112. {
  1113.     char *sine(), *cosine(), *fpdiv(), zero[5];
  1114.     char sresult[5], cresult[5], intres[5], big[5];
  1115.     char *fpmult(), *fpmag();
  1116.  
  1117.     sine(sresult, angle);
  1118.     cosine(cresult, angle);
  1119.     /* check magnitude of denominator :*/
  1120.     /* check magnitude of denominator using exponent */
  1121.     if ((cresult[4] > 128) && (cresult[4] < 132))
  1122.     {
  1123.         initb(big, "30,207,228,127,128");        /*big number */
  1124.         if (sresult[3] > 127)        /*use mantissa sign */
  1125.             return (fpchs(result, big));
  1126.         else
  1127.             return (fpasg(result, big));
  1128.     }
  1129.     /* check for small angle, use small angle approx to
  1130.            avoid underflow */
  1131.     if ((angle[4] < 226) && (angle[4] > 128))
  1132.         return (fpasg(result, angle));
  1133.     else
  1134.         return (fpdiv(result, sresult, cresult));
  1135. }
  1136.  
  1137. /* atanev(result,x) evaluates arctangent for 0 <= x < infinity,
  1138.     result in radians */
  1139.  
  1140. char *atanev(result, x)
  1141.  
  1142. char result[5], x[5];
  1143. {
  1144.     char *coef[8], *piover4, y[5];
  1145.     int index;
  1146.     char *fpadd(), *fpmult(), *atof(), ysq[5], *one;
  1147.     char yterm[5], termreslt[5], *fpasg();
  1148.  
  1149.     piover4 = "\171\356\207\144\0";
  1150.     one = "\0\0\0\100\1";
  1151.  
  1152.     /* small angle approximation */
  1153.     if ((x[4] > 128) && (x[4] < 226))
  1154.         /* use fp exponent to check size, use small angle */
  1155.         return (fpasg(result, x));
  1156.     else
  1157.         fpasg(result, piover4);
  1158.  
  1159.     /* check for argument near one */
  1160.     fpsub(yterm, x, one);
  1161.     /* if it is close to one, as checked by exponent,
  1162.            return the result pi/4  */
  1163.     if ((yterm[4] > 128) && (yterm[4] < 243))
  1164.         return (result);
  1165.  
  1166.  
  1167.     coef[0] = "\270\376\377\177\0";
  1168.     coef[1] = "\277\357\254\252\377";
  1169.     coef[2] = "\106\126\40\146\376";
  1170.     coef[3] = "\146\313\311\270\376";
  1171.     coef[4] = "\34\360\273\142\375";
  1172.     coef[5] = "\206\30\177\215\374";
  1173.     coef[6] = "\154\54\213\131\373";
  1174.     coef[7] = "\67\12\224\275\371";
  1175.  
  1176.     fpadd(termreslt, x, one);
  1177.     fpdiv(y, yterm, termreslt);
  1178.  
  1179.     fpmult(ysq, y, y);
  1180.     /* do poly evaluation */
  1181.     index = 0;
  1182.     /* poly evaluation checked by index limit, and check of
  1183.         variables for under/over flow */
  1184.     while ((index <= 7) && ((y[4] < 100) || (y[4] > 140)))
  1185.     {
  1186.         fpmult(termreslt, coef[index], y);
  1187.         fpadd(result, result, termreslt);
  1188.         index++;
  1189.         fpmult(y, y, ysq);
  1190.     }
  1191.  
  1192.     return (result);
  1193. }
  1194.  
  1195. /* arctan(result,angle) is floating point arctangent evaluation */
  1196.  
  1197. char *arctan(result, x)
  1198. char result[5], x[5];
  1199.  
  1200. {
  1201.     char *atanev(), *fpasg(), y[5];
  1202.     char *fpchs(), *halfpi;
  1203.     int index;
  1204.     halfpi = "\171\356\207\144\1";
  1205.     /* check exponent for very large argument */
  1206.     if ((x[4] > 100) && (x[4] <= 128))
  1207.     {
  1208.         fpasg(result, halfpi);
  1209.     }
  1210.     else        /* go through evaluation */
  1211.     {
  1212.         fpmag(y, x);
  1213.         atanev(result, y);
  1214.     }
  1215.  
  1216.     if (x[3] > 127)
  1217.     {
  1218.         return (fpchs(result, result));
  1219.     }
  1220.     else
  1221.         return (result);
  1222. }
  1223.  
  1224.  
  1225. char *arctan2(result, quadrant, opside, adjside)
  1226.  
  1227. char result[5], opside[5], adjside[5];
  1228. int *quadrant;
  1229.  
  1230. {
  1231.     char x[5], *fpmag(), *fpchs();
  1232.     *arctan();
  1233.     int opsign, adjsign;
  1234.     char *zero, *halfpi;
  1235.  
  1236.     halfpi = "\171\356\207\144\1";
  1237.  
  1238.     zero = "\0\0\0\0\0";
  1239.  
  1240.     opsign = fpcomp(opside, zero);
  1241.     adjsign = fpcomp(adjside, zero);
  1242.  
  1243.     if ((adjsign == 0) && (opsign == 0))
  1244.     {
  1245.         *quadrant = 0;
  1246.         fpasg(result, zero);
  1247.         return (result);
  1248.     }
  1249.  
  1250.     if (((adjside[4] > 128) && (adjside[4] < 226)) || (adjsign == 0))
  1251.         fpasg(result, halfpi);
  1252.     else
  1253.         {
  1254.         fpdiv(x, opside, adjside);
  1255.         fpmag(x, x);
  1256.         arctan(result, x);
  1257.     }
  1258.  
  1259.     if ((adjsign == 0) && (opsign > 0))
  1260.     {
  1261.         *quadrant = 1;
  1262.         return (result);
  1263.     }
  1264.     if ((adjsign == 0) && (opsign < 0))
  1265.     {
  1266.         *quadrant = 4;
  1267.         fpchs(result, result);
  1268.         return (result);
  1269.     }
  1270.     if ((adjsign > 0) && (opsign == 0))
  1271.     {
  1272.         *quadrant = 1;
  1273.         return (result);
  1274.     }
  1275.     if ((adjsign < 0) && (opsign == 0))
  1276.     {
  1277.         *quadrant = 2;
  1278.         fpchs(result, result);
  1279.         return (result);
  1280.     }
  1281.     if ((adjsign > 0) && (opsign > 0))
  1282.     {
  1283.         *quadrant = 1;
  1284.         return (result);
  1285.     }
  1286.     if ((adjsign > 0) && (opsign < 0))
  1287.     {
  1288.         *quadrant = 4;
  1289.         fpchs(result, result);
  1290.         return (result);
  1291.     }
  1292.     if ((adjsign < 0) && (opsign > 0))
  1293.     {
  1294.         *quadrant = 2;
  1295.         fpchs(result, result);
  1296.         return (result);
  1297.     }
  1298.     if ((adjsign < 0) && (opsign < 0))
  1299.     {
  1300.         *quadrant = 3;
  1301.         return (result);
  1302.     }
  1303. }
  1304.  
  1305.  
  1306.     /* page eject */
  1307.  
  1308. /*additional functions for floating point operations*/
  1309.  
  1310. char *sqrt(result, sign, op1)
  1311. char *result, *op1;
  1312. int *sign;
  1313. {
  1314.     char *atof(), *log10(), *fpdiv(), *exp10();
  1315.     char two[5];
  1316.  
  1317.     atof(two, "2");
  1318.     log10(result, sign, op1);
  1319.     fpdiv(result, result, two);
  1320.     exp10(result, result);
  1321.  
  1322.     return (result);
  1323. }
  1324.  
  1325.  
  1326.